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Abstract 

Dark matter (DM) self-interactions have important implications for the formation and evolution of struc- 
ture, from dwarf galaxies to clusters of galaxies. We study the dynamics of self-interacting DM via a light 
mediator, focusing on the quantum resonant regime where the scattering cross section has a non-trivial ve- 
locity dependence. While there are long-standing indications that observations of small scale structure in 
the Universe are not in accord with the predictions of collisionless DM, theoretical study and simulations 
of DM self-interactions have focused on parameter regimes with simple analytic solutions for the scattering 
cross section, with constant or classical velocity (and no angular) dependence. We devise a method that 
allows us to explore the velocity and angular dependence of self-scattering more broadly, in the strongly- 
coupled resonant and classical regimes where many partial modes are necessary for the achieving the result. 
We map out the entire parameter space of DM self-interactions — and implications for structure observa- 
tions — as a function of the coupling and the DM and mediator masses. We derive a new analytic formula 
for describing resonant s-wave scattering. Finally, we show that DM self-interactions can be correlated with 
observations of Sommerfeld enhancements in DM annihilation through indirect detection experiments. 
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L INTRODUCTION 

Dark Matter (DM) is five times as prevalent as ordinary matter, and yet its particle physics 
nature remains elusive. Efforts are underway to detect it through its non-gravitational interaction 
with ordinary matter via direct scattering off nuclei in underground experiments, annihilation to 
Standard Model (SM) by-products in the galaxy today, and direct production in terrestrial col- 
lider experiments. There are well-motivated theoretical reasons to think that DM may reveal itself 
through these means, shedding light on the underlying theory of DM. On the other hand, all ev- 
idence for DM thus far has been obtained through its gravitational interactions, and it remains 
important to continue exploring the nature of DM through its effects on structure in the Universe. 

The formation of structure in the Universe gives critical information about the nature of the 
DM sector. As is well known, the collisionless cold DM (CCDM) paradigm has been highly suc- 
cessful in accounting for large scale structure of the Universe. However, it is far from clear that 
this paradigm can also successfully explain the small scale structure of the Universe. Precision 
observations of dwarf galaxies show DM distributions with cores fl], in contrast to cusps pre- 
dicted by CCDM simulations. It has also been shown that the most massive subhalos in CCDM 
simulations of Miky Way (MW) size halos are too dense to host the observed brightest satellites 
of the MW []2l |3l . Lastly, chemo-dynamic measurements in at least two MW dwarf galaxies show 
that the slopes of the DM density profiles are shallower than predicted by CCDM simulations flU. 

These small scale anomalies, taken at face value, may indicate that other interactions besides 
gravity play a role in structure formation. An interesting possibility is that DM carries self- 
interactions 0. In this scenario, heat can be conducted from the hotter outer to the cooler inner 
parts of the halo through DM collisions, which softens the density profile in the central regions 
of the halo. Recent simulations show that the typical cross section needed to flatten the cores of 
galaxies is a ~ 10~ 24 cm 2 x (mx/GeV) |6]-[8l, where mx is the DM mass. Since it is far larger 
than the typical weak-scale cross section, a ~ 1CT 36 cm 2 , DM candidates cannot be usual WIMPs. 
On the other hand, a light dark force, denoted <p, can provide the required large cross section. A 
perturbative calculation for the scattering cross section gives (in the small velocity limit) 



where the coupling a x is the DM analog of the fine structure constant. Eq. ([T]) shows that light dark 
force with electromagnetic strength coupling can ameliorate discrepancies in small scale structure 
observations. Interestingly, light mediators exist in many DM models which are motivated to solve 
completely different problems [|9]-[T7l. 

Light forces, mediated by a Yukawa potential, can have rich dynamics. The DM self- scattering 
cross section may be velocity dependent, in contrast to the original model where a constant cross 
section is assumed 0. In the regime where axfrix/m^, > 1, Eq. ([T]) breaks down and the non- 
perturbative effect plays a key role in DM scattering. When the momentum transfer is much larger 
than the mediator mass, scattering occurs in the Coulomb limit and the cross section is proportional 
to ~ 1/v 4 with v as the DM relative velocity |fT8l[T9l . While in the quantum resonant regime, the 
scattering cross section can be enhanced and scale as 1/v 2 due to the formation of quasi-bound 
states I1201 . This is the same mechanism that leads to resonant Sommerfeld enhancements in DM 
annihilation [20]. These features have important consequences for DM halo dynamics because 
scattering is enhanced on dwarf galaxy scales compared to MW and cluster scales. It provides a 
natural mechanism to evade constraints from large scales such as elliptical DM halo shapes and 
the Bullet Cluster. It has been shown that most of parameter space of interest for thermal DM is in 
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this non-perturbative regime [|20ll . 

While self-interacting DM has been the subject of astrophysical interest and numerical simu- 
lation, the particle physics aspects have been comparatively little examined. Studies have so far 
limited themselves to regions that can be approximated analytically through classical or Born for- 
mulae ll2Tl - |25l , or else have considered a limited range of parameter space [|26ll . The purpose of 
this paper is to study the full range of effects of light dark force dynamics on halo structure, and 
is intended as a companion paper to ll20ll . In ll20ll . we have laid out a model of self-interacting 
DM that satisfies relic density considerations while giving rise to a rich structure in the scattering, 
including the presence of velocity dependent resonances. Here we delve into many more details. 
We discuss the method that we use for an improvement in the numerical efficiency for solving the 
Schrodinger equation such that we are able to reach regions of scattering parameter space where 
many partial wave i modes are required. This method allows us to explore the strongly-coupled 
resonant and classical regimes. We are able to verify numerically the classical formula, which has 
never been done before. We are also able to examine the angular dependence of the scattering 
cross section in the classical and strongly-coupled regimes, observing the transition to the weakly 
coupled regime with forward-peaked Rutherford scattering. We examine in detail s- and p-wave 
resonances in the strongly-coupled regime, and provide benchmark points for simulations. 

The outline of this paper is as follows. First, we discuss the case for self-interacting DM, sum- 
marizing the current status of DM simulations and observations of small scale structure. Then we 
describe our setup, diving into technical details of solving the dynamics of strongly-coupled sys- 
tems. We show the results of our method, encapsulating the Born, resonant, and classical regimes, 
and we examine the velocity and angular dependent scattering effects on halo structure. We then 
present a new analytic result for the s-wave resonant regime before connecting our method to relic 
density calculations, explicitly including the effect of the Sommerfeld enhancement. Lastly, we 
discuss connections to observation in indirect detection experiments, showing how the enhance- 
ment in self-scattering can also be important for DM annihilation. We then conclude. 

H. SELF-INTERACTING DARK MATTER AND SMALL SCALE STRUCTURE 

For some time there has been debate about whether the paradigm of CCDM accurately de- 
scribes the observed small scale structure in the Universe. Small scale objects (e.g., dwarf galax- 
ies) are typically DM-dominated, and therefore offer potentially cleaner laboratories to test CCDM 
predictions compared to systems with higher baryon densities. Here, we describe three discrepan- 
cies, and show how two of them may point beyond CCDM in the form of DM self-interactions. 
We emphasize, however, that the situation remains far from clear, and ultimately more detailed 
numerical simulations including baryonic effects are required before drawing definitive conclu- 
sions II271I28TI. 

Core-vs-cusp problem: The central density profiles of dwarf galaxy halos indicate a long- 
standing discrepancy between steep cusps predicted by CCDM-only simulations Il29l - |3"TI com- 
pared to flat cores inferred from observed galaxy rotation curves [Q] l32l - [35ll . Observations of 
clusters of galaxies may also exhibit cored profiles ll3~oT - [38ll . Baryonic effects may provide an as- 
trophysical mechanism for flattening the DM density profile in the center of a galaxy (or cluster of 
galaxies), which are often baryon dominated. It has been argued that feedback from (dissipational) 
baryonic matter leads to further contraction of the central DM cusp [391, further exacerbating the 
discrepancy. However, simulations have shown this mechanism to be less effective than previ- 
ously thought BOl |4TT| . Moreover, supernova feedback may have the opposite effect: supernova 
energy injected into the interstellar medium leads to baryonic outflow, which can gravitationally 



3 



disrupt the central cusp, resulting in lower DM densities compared to CCDM-only simulations 
Il42l - l47ll . However, this mechanism seems unlikely to explain central cores in metal-poor galaxies 
with limited star formation rates [[341 . 

Missing satellites problem: There has been an order of magnitude discrepancy between the 
number of observed and expected satellites of the Milky Way (MW) ||48] - [5T|I . Baryonic processes 
such as supernova feedback and/or photoionization may play important for suppressing star forma- 
tion in dwarf galaxies, explaining the observed (weak) baryon content of these small galaxies ll52l . 
Recently, the Sloan Digital Sky Survey has discovered many faint galaxies, such that it is evident 
that as many as a factor of 5 — 20 of the known dwarf galaxies could be still undiscovered due to 
faintness, luminosity bias and limited sky coverage Il53l - l55l . Consensus is thus shifting toward the 
view that the number of MW subhalos is not an issue for the predictions of CCDM, at least for 
smaller subhalos. 

Too-big-to-fail problem: Detailed studies of the brightest MW dwarf spheroidal (dSph) galax- 
ies, which are DM dominated at all radii, show discrepancies with CCDM-only predictions (see 
e.g. [56J). These satellites are expected to be hosted by the largest subhalos in the MW halo, as 
they have the largest velocity dispersions observed from their rotation curves. However, the most 
massive subhalos predicted by CCDM-only simulations are too massive, with central densities too 
large, to host the brightest observed satellites [5Vj . Simulations predict O(10) subhalos with 
maximum circular velocity V max > 30 km/s, whereas the MW dSphs have V max < 25 km/s. 1 This 
discrepancy may share a common resolution with the core-vs-cusp problem; these massive subha- 
los can be reconciled with the observed dSphs if their central densities are reduced compared to 
CCDM predictions. Indeed, analyses of stellar subpopulations within several dSphs indicate cored 
central density profiles Iffl I58W6TI (except for Draco [62]). Several studies using hydrodynamical 
simulations have suggested that baryonic physics — i.e., feedback from star formation and su- 
pernovae, as well as ram pressure and tidal stripping from the host halo — may induce dSph 
cores ll63l - [66ll . while Ref. 11671 found a smaller impact from baryonic effects. Additionally, the 
severity of this problem can be reduced by taking into account statistical variation in the formation 
of MW-sized halos [68], as well as uncertainty in the MW halo mass which sets the normaliza- 
tion of the subhalo mass spectrum. Larger MW halo masses lead to a larger discrepancy between 
simulation and observed MW satellites. For example, Ref. 11641 used a MW mass of 8 x 1O U M 
and saw no too-big-to-fail problem. On the other hand, Ref. ||69ll argued that a larger MW mass, 
around 2 x 1O 12 M , is warranted. Whether this larger estimate or the smaller one advocated in 
[TTOl prevails will have important implications for the too-big-to-fail problem. 

Given these persistent questions about the accordance of observations with the predictions of 
CCDM, it is interesting to look beyond the paradigm of cold and collisionless DM. One of the first 
attempts to do this was to give DM some kinetic energy, i.e., to make it warm. Warm DM predicts 
a suppression in the halo mass function at small scales, below the free- streaming length. Thus, 
warm DM effectively removes substructure, and predicts a reduced number of satellites in a galaxy 
such as the MW. On the other hand, warm DM halos may be less concentrated than CCDM halos 
on scales of order the free- streaming length, but they are still cuspy iTTTl l72l . As a result, warm 
DM solves only the missing satellites problem, which is considered the least severe discrepancy, 
but not the remaining problems. 

The other known mechanism for changing the structure of DM halos is self-interactions. Self- 
interacting DM was introduced as a solution to the core-vs-cusp and missing satellites problems 



These predicted most massive subhalos are "too big to fail" in forming stars, unlike shallower potentials of much 
smaller subhalos. 
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in Ref. 0. Self-interactions cause energy transfer from the hotter outer halo to the colder central 
region, thereby forming a core. At the same time, collisional stripping of dwarf subhalos within 
the hotter MW host halo can deplete the abundances of satellites. Early simulations, which fo- 
cused primarily on the case of a constant (velocity-independent) scattering cross section, found 
that a /m x ~ 1 — 10 cm 2 /g flattened the central densities in dwarf galaxies in accordance with 
observations and cr/mx ~ 10 cm 2 /g reduced significantly the number of MW subhalos 11731 . 

Subsequent studies, however, found rather serious problems with self-interacting DM due to 
conflicts with other observations. The simulation of Ref. [74] concluded that cr/mx < 0.1 cm 2 /g 
is required to avoid core formation in cluster halos in conflict with gravitational lensing obser- 
vations of CL 0024+1654. Ref. 11751 argued that o /m x < 0.02 cm 2 /g is required by cluster 
ellipticity constraints, while Ref. 11761 showed that a /nix ~ 0.3 — 10 4 cm 2 /g is excluded by re- 
quiring that elliptical galaxy halos do not evaporate within hot cluster halos. Lastly, Ref. [17711 
obtained a/m x < 1 cm 2 /g from the X-ray and lensing observations of the Bullet cluster. 

More recently, there have been two major developments leading to a revival of self-interacting 
DM. First, DM self-interactions need not have a cross section that is constant in velocity lfT8l - l2~3l . 
For light dark force mediators, once the momentum transfer becomes comparable to the mediator 
mass, the cross section begins to decrease rapidly (analogous to Rutherford scattering). Since 
larger halos have larger characteristic velocities, the cross section can be large in dwarf galaxies 
(v ~ 10 km/s), but negligible on cluster scales (v ~ 1000 km/s) to evade the aforementioned 
constraints. 

Second, considerable progress has been made in numerical simulations of self-interacting 
DM ||6T]8ll78l. In particular, the issue of self-interacting DM constraints from galaxy clusters was 
recently revisited in Refs. B71|78l. In these simulations, a very different conclusion was reached 
from earlier simulations. In particular, the constraints from cluster halo triaxiality were found to 
be much weaker than previously estimated. They conclude that previous works did not take into 
account that the observed ellipticity has contributions from regions well outside the core, and this 
region retains its triaxiality. They also find that the residual triaxiality is larger than previously 
estimated 11731 . and that the remaining discrepancy can be accounted for in the ellipticity scatter 
between different DM halos. Furthermore, the authors also find that the tendency of subhalos 
to evaporate is not significant for a/m ~ 1 cm 2 /g. Lastly, the cluster CK 0024+1654 used by 
Ref. [74] is now known to be undergoing a merger along the line of sight, making it less useful as 
a comparison case with non-merging simulation data. 

Overall, while the situation for self-interacting DM is not yet resolved, much progress has been 
made. The most recent simulations have shown that a/mx ~ 0.1 — 10 cm 2 /g on dwarf scales 
is sufficient 2 to solve the core-vs-cusp and too-big-to-fail problems [[6T-l8l 1781. while constraints 
on MW and cluster scales require a/m x < 0.1 — 1 cm 2 /g [[6l [71 1751 . It appears that all the 
data may be accounted for with a constant scattering cross section around cr/m x ~ 0.5 cm 2 /g. 
On the other hand, particle physics models of self-interacting DM generically predict a velocity- 
dependent scattering cross section over a wide range of parameter space, as we discuss below. 



2 Ref. [8| found that a/mx = 0.1 cm 2 /g is too small, although the precise lower bound is unknown. 
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III. DARK FORCES AND DARK MATTER SCATTERING 



In order to explain astrophysical observations on dwarf galaxy scales, the DM elastic scattering 
cross section must be 

a ~ 1 cm 2 (m x /g) « 2 x 1(T 24 cm 2 (m x /GeV) , (2) 

which is much larger than a typical weak-scale cross section a ~ 10 -36 cm 2 . Therefore, this 
suggests the existence of a dark force boson <p that is much lighter than the weak scale. 

In this work, we consider a phenomenological approach where nonrelativistic DM scattering is 
described by a Yukawa potential 

V(r) = ±— e- m f\ (3) 
r 

which can be either repulsive (+) or attractive (— ). This interaction arises for as a vector or 
scalar mediator, with interaction 



^ ( g x X'j^X(j) fl vector mediator 
^ gxXX<p scalar mediator 



and dark fine structure constant ax = 5 , x/(4vr). Scalar interactions are purely attractive, while 
a vector interaction is both attractive (XX scattering) and repulsive (XX or XX scattering). 
Thus, in the vector case, asymmetric DM (X only) will have purely repulsive interactions, while 
symmetric DM (equal X, X) will have both attractive and repulsive interactions, with the total 
effective cross section given by the average of the two. 

Numerical N-body simulations have investigated the impact of DM self-interactions on struc- 
ture formation. The relevant input is the differential cross section da/dVt, as a function of the 
DM relative velocity v. Since simulations track particle trajectories before and after collisions, 
the angular distribution over the scattering angle 9 is important. However, to compare across 
different parameter regions, with different angular dependencies, it is useful to consider an inte- 
grated cross section that captures the relevant physics. The usual quantity is the standard cross 
section a = J dVt(do / dVt) . However, for light mediators, a receives a strong enhancement in the 
forward-scattering limit (cos 9 — > 1), and for the purposes of affecting the DM distribution this 
enhancement is spurious since the DM particle trajectories are unchanged. In the plasma literature, 
two additional cross sections are defined to parametrize transport IT791I , the transfer cross section 
a t and the viscosity (or conductivity) cross section ay: 

dQ (1 - cos 9) , ay = / dQ sin 2 9-^-. (5) 

The transfer cross section is weighted by (1 — cos 9), the fractional longitudinal momentum trans- 
fer, while the viscosity cross section is weighted by the energy transfer in the transverse direction, 
sin 2 9. The transfer cross section has been used in the DM literature to regulate the forward- 
scattering divergence. On the other hand, the viscosity cross section weighs forward and back- 
ward scattering evenly. It takes into account that forward and backward scattering affect the DM 
halo equally, since DM particles simply exchange trajectories that they would have had in the ab- 
sence of a collision. It also takes into account that we expect that perpendicular scattering is most 
efficient for "thermalizing" the DM halo and affecting structure observables. 

In addition, the transfer cross section obviously fails if DM scattering occurs between identical 
particles. Taking quantum indistinguishability into account, both forward and backward scattering 
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diverges, corresponding to poles in the t- and w-channel diagrams, ctt regulates only the for- 
ward divergence, making it an inadequate description for the case of quantum indistinguishable 
particles. Since both forward and backward scattering leave the DM distribution unchanged, the 
relevant cross section should regulate both divergences, which ay does, but ot does not. 

In order to make contact with previous work, however, we focus on a T , rather than ay. Un- 
der the assumption of classical distinguishibility in scattering, we find that ot and ay differ by 
less than a factor of two, with ay for distinguishable and indistinguishable particles differing by 
another 0(1) number. Thus the overall effect both of distinguishability and of the transfer versus 
viscosity cross section is (9(1). For the purpose of presenting our results, we assume classical dis- 
tinguishability and take a T as a suitable measure for the effects of DM scattering on halo shapes. 
Of course, a full-scale N-body simulation should make use of the angular information in the dif- 
ferential scattering cross section, da /dVL, and do away with the proxy of a transfer or viscosity 
cross section altogether, though in most cases the difference between the results using ay or a T 



versus da /dVt will be small. In Sec. IV C below, we discuss the angular dependence in more detail 
and present benchmarks for simulation. 

The transfer cross section, computed perturbatively in ax from Eq. (|3]), is given by 



On 



= log 1 + m\v 2 ml) - x 2 - . (6) 

m x v 4 V v mi+m x v 2 J 

for both attractive and repulsive potentials lETI . where v is the relative velocity. This perturbative 
expression is valid only within the Born approximation, requiring axmx/m^ <C 1. Outside 
this limit, the Born approximation is not valid and non-perturbative corrections become crucially 
important. 

Within the non-perturbative regime, analytic formulae for a T have been obtained only within 
the classical limit (mxv/m^ 3> 1) ll2Tll8~0ll8T1l . given for an attractive potential by 



_clas 

a T 



l^ln^ + r 1 ) /? < 10- 1 

f!/3 2 / (1 + 1.5/3 1 - 65 ) 10- 1 < /? < 10 3 (7) 



(ln/3 + 1 - iln- 1 /?) 2 /3>10 3 



where /3 = 2ax'm ( f ) /(mxv 2 ). Many previous works Il2fl - l23ll25l . including recent N-body simula- 
tions [6], have focused specifically on the case where DM scattering is described by an attractive, 
classical cross section, given by Eq. ([7]). We emphasize, however, that this case is just one out 
of many possibilities, and in general the non-perturbative regime remains largely unexplored. We 
collect, for reference, the analytic formulae in the Appendix for the Born, attractive and repulsive 
classical, and s-wave resonance cases. 

For a large parametric range of interest for DM self-interactions, both quantum mechanical 
and non-perturbative effects become important, and neither the Born nor classical approxima- 
tions are valid. The onset of these effects is governed by the conditions axTnx/m^ > 1 and 
mxv/m^ < 1, respectively. We denote this region of parameter space as the "resonant regime," 
since one important effect is the appearance of quantum mechanical resonances in ax correspond- 
ing to (quasi-)bound states in the potential. 

Within the resonant regime, there exists no analytic formula for a T , and it must be computed by 
solving the Schrodinger equation directly using a partial wave analysis. The differential scattering 
cross section is given by 



da 

dn k 2 



-I ^2(2£ + l)e i5t P e (cos9) sin 5* 



(8) 



1=0 
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where 5? is the phase shift for a partial wave i. In terms of the phase shifts, the transfer cross 
section is given by 

ark 2 



Air 



£^+l)sin 2 (fc + i-fc). (9) 

1=0 

To obtain Sg, one must solve the Schrodinger equation for the radial wavefunction Ri(r) for the 
reduced DM two-particle system, given by 

1 d / 2 dRi\ r,-, i(£+l 



r 2 dr V dr 



( k2 ~ ~ ^V(r))Re = (10) 



with reduced mass /i = m x /2 and momentum k — fxv. The phase shift 5i parametrizes the 
asymptotic solution for Ri(r), given by 

lim R(.{r) oc cos5eje(kr) — sin. St ne(kr) , (11) 

r— >oo 

where jg (n e ) is the spherical Bessel (Neumann) function. 



IV. NUMERICAL SCATTERING RESULTS 

In this section, we present our numerical results. First, we describe our numerical method for 
computing the DM self-interaction cross section or. Next, we investigate the velocity-dependence 
and angular-dependence of DM scattering. For realistic particle physics models of self-interacting 
DM, scattering can possess a wide range of nontrivial dependence on velocity and scattering an- 
gle, whereas N-body simulations have considered isotropic scattering with constant or particular 
choices of velocity dependencies. 



A. Numerical Method 

To solve the Schrodinger equation, it is useful to define the variables ll26l 

Xe = rR e , x = a x m x r, a=- — , b= , (12) 

2a x 



such that Eq. ( 10 ) can be expressed as 



Ml + a2 _^+i) ± i e _, /> h<(l) = . (13) 



\ dx 2 x 2 x 

To compute or, we first compute 5e for given (a, b, €) as follows. 

1. We impose an initial condition for xi an d x'i at a point x = Xi close to the origin. For 



Xi <C b, (£ + 1)/ a, Eq. ( fl~3] ) is dominated by the angular momentum term, and we expect 
Xi( x ) x £+1 . Thus, we take Xe( x i) — 1 an d Xe( x i) = + l)/ x i\ the overall normalization 
is irrelevant. 



2. We solve Eq. ( [13] ) numerically within the domain Xi < x < x m . The matching point x m is 
determined by the condition a 2 exp(— x m /b)/x m , where the potential term is suppressed 
compared to the kinetic term. 
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b = a x m x /m^ b = a x m x lm (P 



FIG. 1: Colored regions show parameter points (a, b) within our numerical scan, with the corresponding 
values of ark 2 / (Air) (left) and £ max (right) at each point. The classical, Born, and resonant regimes are 
delineated by solid lines. 



3. At match xe ( an d its first derivative) onto the asymptotic solution, given by 

Xt(%) oc xe l&l (cos^ je(ax) — sin^n^(ax)) . (14) 
Inverting Eq. (TBI), the phase shift is given by 



(ll ffl j|((ll m ) (3(> j £ (ciX m) ^mXii-^m) 1 /ic\ 
tari()£= r r , P e = r 1 (15) 

ax m n' t (ax m ) - p t n e (ax m ) Xe{ x m) 

in terms of our numerical solution for xe at x m- Our numerical method makes an initial 
guess for (xi,x m ) and computes 5e, and then successively decreases (increases) Xi {x m ) 
until 5e converges at 1%. 



max 



4. The last step is computing by summing Eq. (|9|) over £, truncating at £ max . We iterate 
until ot converges to 1% and 5^ max < 0.01 through ten successive iterations. This condition 
is quite conservative, typically summing many more £-modes than required. 

For a given (a, b), we can then express ot in terms of the physical parameters (m x , m^,, ax, v). 
Our numerical code for this solution was written using Mat hematic a. 

With our numerical method in hand, we performed a fine-grained scan over 2 x 10 5 parameter 
points (a, b). Fig. [T] gives a birds-eye view of our full numerical dataset, with the colored points 
showing the parameters (a, b) in our scan. In the left panel, the different colors correspond to the 
computed value of ark 2 / (4tt) obtained from Eq. ([9]), with the corresponding value of £ max shown 
in the right panel. The white region (upper right) was omitted from our scan. The solid lines at 
b = 1 and 2ab = 1 delineate the Born regime (b 1), the classical regime (2ab ^1), and the 
resonant regime (b > 1 and 2ab < 1). The latter exhibits a pattern of resonances in cr^. 

There is an importance difference between our method and that of Ref. 11261 . which performed 
a similar calculation of a T within the resonant regime, albeit for a limited choice of parameters. 
Ref. [26J obtained Si by matching onto an asymptotic form Xi{ x ) ^ sin(ax — tt£/2 + 5i), which 



is equivalent to Eq. ( |T4| ) for x — > oo. For finite x, this form is valid only if both the Yukawa 



and angular momentum terms in Eq. ( 13 ) are suppressed compared to the kinetic term, whereas 



Eq. ( 14 ) requires only the Yukawa term to be suppressed. Therefore, as i is increased, the method 
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FIG. 2: Left: Numerical calculation of ar/mx, truncated at fixed l max , showing convergence with in- 
creasing ^ max - The parameter point chosen corresponds to the classical regime with an attractive potential. 
The convergence to the classical analytic result shown by dashed line. Right: Numerical calculation (solid 
blue) of &t I m x versus m^, showing convergence to the classical analytical formula (dotted pink) and Born 
approximation (dashed gold) in the classical and Born regimes. 



of Ref. [26J requires integrating Eq. ( [13] ) to much larger x than in our method, and is therefore 
much less efficient. Thus, Ref. [26] truncates at £ max = 5 in their calculation, whereas we are 
able to perform efficient calculations with £ max ~ 1000. We demonstrate this point in Fig. [2[ 
showing how ot depends on £ max for one parameter choice in the classical regime. Our numerical 
calculation (solid line) converges for £ max > 1000, in good agreement with the classical cross 
section (dashed line). 3 

We can also see the convergence to classical and Born analytic formulae in the right panel of 
Fig.[2j The dashed gold and dotted pink lines show the results for the Born and classical analytic 
formulae, and we see that in the regime of validity, our numerical results (solid blue line) agree 
well with the analytic formulae. In the quantum resonant regime, neither of the analytic formulae 
reproduce the behavior of the resonant peaks and anti-resonant valleys. Also note that the Born 
approximation over-estimates the cross section in the classical regime. 



B. Velocity-dependence in dark matter scattering 

The most important feature that emerges from our numerical study is the highly nontrivial 
velocity-dependence of o T within the resonant regime. While previous studies have focused on 
either constant ot or specific v -dependencies, a rich array of possibilities can arise in general, and 
the velocity behavior can be rather complicated. 

In Fig. [3] we show the cross section as a function of velocity for an attractive potential with 
ax = 10 . Each curve corresponds to a different value for b (where b = ax^x/rn^), as 
indicated by the numerical values in the figures. The quantity a T m\ is a useful normalization 
for the cross section since, for fixed ax, it depends on v and mx/m^ only (as opposed to mx 
and separately). Thus, to obtain the required level scattering in dwarf halos, each curve can 



3 The reader should not be troubled by the fact that a? can be negative for certain values of i max . Due to the fact 
that the momentum and orbital angular momentum operators do not commute, the transfer cross section, defined in 
terms of momentum eigenstates, is a physical quantity only in the limit £ max — > oo, not for a particular value of t. 
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FIG. 3: Velocity-dependence of the scattering cross section for an attractive potential with ax = 10 , 
computed numerically (solid lines), for various values of b labeling each curve. Numerical values indicate 
axmx /mj, chosen for each curve. Dashed lines show extrapolation using classical formulae. Each curve 
can be normalized to ax/mx ~ 1 cm 2 /g on dwarf halo scales (v « 10 km/s) by dividing by m x for a 
particular choice of mx- 

be normalized to a T /m x ~ 1 cm 2 /g at v m 10 km/s by choosing m x appropriately, which also 
fixes m^. 

The cross sections shown in Fig. [3] exhibit a wide variety of behaviors and features. The se- 
quence of different cross sections, ascending from b = 1 to b ~ 1000, illustrate the onset of 
resonance features beyond the Born regime (b <C 1). Increasing b, we first see the appearance of 
an s-wave resonance for b = 1.68; the phase shift behaves as |d>o| — > tt/2 for v — > 0, and so the 
cross section becomes strongly enhanced, growing as o T — > 16% / (m x v 2 ) on-resonance. Moving 
to larger values of b, the cross section becomes reduced, and we see the appearance of an s-wave 
antiresonance for b = 4.52, where the cross section is strongly suppressed at low velocity. Higher 
£-mode resonances appear as peaks at finite v where o T is enhanced. For b = 17.6 we note the 
appearance of a p-wave resonance at v ~ 30 km/s, and for higher b, spectral features become 
increasingly prevalent. At high velocity, all cross sections converge to the same Coulomb result, 
gtvp^x v ~ 4 i independent of m^/mx- 

In Fig. |4} we show a similar set of results for the cross section arising from a repulsive inter- 
action, with ax = 10~ 2 . Unlike the attractive case, no resonances arise in the "resonant" regime, 
since there are no bound states in the potential. However, the cross section exhibits a clear velocity 
dependence where scattering on dwarf scales can be enhanced compared to larger scales. Larger 
values of b (i.e., smaller m^, for fixed m x ) correspond to a longer range force, enhancing ax- 

C. Angular dependence in dark matter scattering 

Since N-body simulations track particle trajectories before and after collisions, it is required to 
know the differential cross section and its dependence on the scattering angle 9. Although for s- 



4 Although the cross sections shown in Figs.|3]and|4]have fixed ax = 10~ 2 , these results can be generalized to other 
values of ax by a shift in the horizontal and vertical axes. Effectively, this shift amounts to relabeling the axes by 
v — > v x (10~ 2 /a x ) cinda T m x — > a T rr? x x (a x /10~ 2 ) 2 . 
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wave scattering, da/dVL is isotropic, more complicated angular dependencies arise in a wide range 
of parameter space. In general, if DM scattering is velocity-dependent, then often the differential 
cross section carries a nontrivial angular dependence. 

First, we investigate the impact of anisotropic scattering within the classical regime. The nu- 
merical simulation of Ref. [6J considered a velocity-dependent cross section given by Eq. (Q, 
corresponding to an attractive interaction. Here, we consider one specific benchmark point from 
this work (denoted therein as "RefP2"), shown to solve small scale structure anomalies. This 
benchmark is parametrized phenomenologically by cr™ ax /mx = 3.5 cm 2 /g and t> max = 30 km/s; 
these quantities are related to the underlying parameters (m x , m^, ax) by a™ 3 * = 22.7 /m 2 ± and 
^max = ^ m 4> a x/(^rnx)- 5 We emphasize that Ref. [|6] assumed in their simulation an isotropic 
differential cross section given by da/dVL = o"r/(47r). With our numerical solution in hand, we 
can check whether this approximation is justified. 

In Fig. [5] (left), we present our results for da /dVL for the RefP2 benchmark point, with each 
panel corresponding to a different velocity. The horizontal black lines show the isotropic ap- 
proximation da/dVL = a T /(An) adopted by Ref. |5|. The solid curves show our numerical cal- 
culation of da/dVL. Although t> max and cr™ ax /mx are fixed, an additional input is required to 
fix the three parameters (mx,m^, ax)- We have taken mxv/m^ = 10 (thick blue curve) and 
mxv/rricf, = 50 (thin green curve) 6 ; to the extent that these curves overlap, da/dVL does not 
depend on this additional parametric freedom. The dashed red line shows the usual Rutherford 
formula da /dVt = a\j{w? x vi sm4 #/2). From these plots, we conclude: 

• At small velocity, da/dVt has a nontrivial angular dependence, with many small-scale angu- 
lar features oscillating about a nearly fiat profile. Since astrophysical structure observables 
are likely insensitive to small-angle features, we conclude that the isotropic approximation 
appears well-justified in this regime. This behavior is expected since (3 3> 1 corresponds to 
the strong coupling limit, and the Yukawa potential approaches the hard sphere limit with 
radius set by mT , with da /dVt fiat. 

• At large velocity, da/dVt becomes peaked for forward- scattering (cos 6 — > 1). This behavior 



5 To clarify this notation, we note that the quantity a^v is maximized at v — u max , at which cry = cr™ ax ||23l . 

6 For visual clarity, we have smoothed these curves by averaging each point over an interval Acos# — ±0.01 to 
eliminate small-angle features. 
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FIG. 5: Left: numerical solution for da/dQ m^ 1 for Ref. [6] benchmark point for mxv/m^ = 10 (thick 
blue) and 50 (thin green). Right: numerical solution for da /dSlrrf^ for benchmark point with p-wave 
resonance atw rj 10 km/s (solid blue). Exact results are compared to the isotropic approximation da/dQ = 
<7r/(47r) (flat black) and the Rutherford formula (dashed red). 
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is expected since /? <C 1 corresponds to the Coulomb limit, and da / dfl is well-approximated 
by the Rutherford formula. We conclude that an isotropic approximation is not valid in this 
limit. However, since the cross section is suppressed at larger velocity, this discrepency may 
not be important. 

Similar conclusions apply to other parameter points in the classical regime: scattering is approxi- 
mately isotropic for /3 > 1, but becomes forward-peaked for (3 < 1. 

Next, we consider a benchmark parameter point within the resonant regime: mx = 100 GeV, 
m<f> = 17 MeV, ax = 3 x 1CT 3 . These parameters have been chosen to give a p-wave resonance on 
dwarf scales, with a peak at v — 10 km/s with ax/mx = 22.5 cm 2 /g. In Fig. [5] (right), we show 
our numerical results for do /dVt (solid blue curves), with each panel corresponding to a different 
velocity, compared to the isotropic approximation da/dVL = 0" T /(47r) (horizontal black lines). 
At small v, scattering is predominantly s-wave, with da/dVL nearly flat. At v — 10 km/s, the 
£ = 1 term dominates, enhancing the scattering cross section and giving an angular dependence of 
da/dVL oc cos 2 9. For larger v, higher £ modes become important, and da /dVt becomes forward- 
peaked, approaching the Rutherford formula. For a p-wave resonance, it is clear that the angular 
dependence is crucial. Although ar/mx may be strongly enhanced on dwarf scales, the impact 
on astrophysical structure observables is likely less pronounced. The p-wave angular distribution 
is weighted toward forward or backward scattering, whereas we expect structure observables to be 
more sensitive to perpendicular scattering (cos 9 ps 0). 



V. PARAMETER SPACE FOR SELF-INTERACTING DARK MATTER 

In this section, we show how bounds from astrophysical observations of structure map onto 
the underlying DM particle physics parameter space. Within our simple framework, there are 
only three parameters (mjf,m^ as well as one overall sign corresponding to a repulsive or 
attractive potential in Eq. Q. For a given parameter choice, we compute (Jt(v) either numerically 
or using the Born or classical analytic approximations (where valid). However, the DM scattering 
probability within a halo is determined not by one fixed v, but rather by a convolution over different 
velocities and densities as a DM particle traverses the halo, requiring detailed N-body simulations 
which are beyond the scope of our work. Instead, we consider the velocity-averaged transfer 
cross section (a T ) as a suitable proxy for the quantity being constrained by astrophysical bounds. 
Averaging over the initial DM velocities v 12 with exponential weight, we have 

M - / Mlft " 41) - / e '^'° ° T{V) ' (16) 

where t> is the most probable velocity and v = \vi — is the relative velocity. We choose vo to 
be characteristic of different size halos, described below. Although velocity-averaging is clearly 
irrelevant for a constant cross section, it is especially important for strongly velocity-dependent 
cross sections (e.g., resonant features). 

Our results for (cr T ) are presented in Fig. [6j For both attractive (left) and repulsive (right) 
potentials, we show the allowed range of (mx,m^) for ax = 10 _1 , 10~ 2 , 1CT 3 . Astrophysical 
bounds on different scales are indicated as follows: 

• Blue regions show 0.1 < (<y T /m x ) < 1 cm 2 /g (light) and 1 < (a T ) < 10 cm 2 /g (dark) 
on dwarf scales (t>o = 10 km/s), required for solving small scale structure anomalies. 

• Red contours show {o T )/m x = 0.1 and 1 cm 2 /g on MW scales (v = 200 km/s). 
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• Green contours show (ax) /mx = 0.1 and 1 cm 2 /g on cluster scales (vo = 1000 km/s). 

The dashed lines indicate where we use analytic formulae for a T , given in Eq. ([7]), to interpolate 
our results into the classical (top) and Born (bottom) regimes. The fact that our numerically com- 
puted contours match well onto these regimes demonstrates the consistency between the numerical 
and analytic results. 

Since N-body simulations have been performed for only a limited choice of cross sections, 
the precise numerical values of these constraints are open to interpretation. For a constant cross 
section, Ref. [7 J found that ax/mx = 1 cm 2 /g is too large, causing too-small central densities 
in dwarf spheroidals and clusters and is marginally consistent with ellipticity constraints on MW 
scales, while ax/m x = 0.1 cm 2 /g satisfies all constraints including on dwarf scales. On the 
other hand, simulations with a velocity-dependent cross section (assuming a classical, attractive 
form for ax) have favored larger values on dwarf scales, a T /m x ~ 10 cm 2 /g @. Therefore, 
we expect that the actual astrophysical bound on MW (cluster) scales lies between the red (green) 
lines between (crx)/mx = 0.1 — 1 cm 2 /g, with the area to the left excluded, while the blue 
regions are favored by solving small scale structure anomalies. More precise limits require future 
N-body simulations utilizing the full velocity-dependent form for ax(v), as a function of the DM 
parameters. 

The most striking features emerging from our numerical calculation are the pattern of quantum 
mechanical resonances and antiresonances for the attractive potential case (absent for the repulsive 
case). For fixed (cr T ) /m x , the (anti)resonances favor larger (smaller) m x , corresponding to peaks 
pointing to the upper right (lower left) in Fig. [6j These features are more pronounced for smaller 
v and larger ax since the conditions rrixv/m^ < 1 and ax m x / m <p ?t 1 govern the onset of 
quantum mechanical and non-perturbative effects, respectively. It is clear that the resonant regime 
corresponds to a large region of parameter space, mx ~ GeV — TeV, where ax is computed 
numerically. In the next section, we will derive an analytical formula for ax in the resonant 
regime. 

Our main conclusion from Fig. [6] is that for a wide range of (ax,m x ), self-interacting DM 
can explain small scale structure anomalies while remaining consistent with other astrophysical 
bounds. 

• A wide range for the DM mass m x is allowed, from sub-GeV to multi-TeV or beyond. 

• A wide range of perturbative couplings ax are allowed; we explicitly showed results for ax 
between 10 _1 and 10~ 3 . 

• For fixed (mx, ax), the mediator mass is determined within an order of magnitude by as- 
trophysical bounds. Generally, for mx < TeV, we require m^ ~ 1 — 100 MeV, with smaller 

777.0 for mx > TeV. 

Future observations on MW and cluster scales can play a key role in narrowing this parame- 
ter space by giving additional velocity data points for ax(v). For example, evidence for self- 
interactions on larger scales at the level of (ax) /mx ~ 0.1 cm 2 /g would favor light DM at the 
GeV-scale; excluding self-interactions below this level would favor heavier DM. 

VI. RESONANT s-WAVE SCATTERING: ANALYTIC RESULTS 

We derive a new analytic formula for the s-wave scattering cross section that is valid in the 
resonant regime. This result provides an accurate description of DM scattering in a parameter 
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FIG. 6: Parameter space consistent with astrophysical bounds for attractive (left) and repulsive (right) poten- 
tials for different ax- Blue regions show where DM self-scattering solves small scale structure anomalies, 
while red (green) show bounds on Milky Way (cluster) scales. Numerical values give (ax) /mx in cm 2 /g 
on dwarf ("dw"), Milky Way ("MW"), and cluster ("cl") scales. See text for details. 
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region which has not been previously analytically accessible and is complementary to the Born 
and classical regimes. We give a simple analytic condition for resonances and antiresonces to 
occur, and we confirm our results against our numerical computation. 

Although the Schrodinger equation cannot be solved analytically for the Yukawa potential in 
the non-perturbative regime, a useful proxy is provided by the Hulthen potential 



a x 5 e~ 5r 
1 -e" 



which is analytically solvable for £ = 0. The Yukawa and Hulthen potentials behave similarly, 
scaling as 1/r at short distances and becoming screened for large distances. The Hulthen screening 
mass 5 is assumed to be related to by 5 = nm^, where k is an 0(1) numerical constant. In 
computing the Sommerfeld enhancement for DM annihilation, Ref. ||871 showed that Eq. ( fT7 ) 
provides an accurate analytic approximation of the Yukawa potential, with k = n 2 /6. Here, we 
follow a similar analysis to compute the cross section for DM scattering; however, we keep k as a 
free parameter. 

Defining c = axmx/5 and substituting the Hulthen for Yukawa potentials, Eq. ( |T~3| ) becomes 

d 2 



&+°'*T^)x>W = > (18) 

for i = 0. With another change of variables t — 1 — c~ x l c and Xo( x ) = tO- ~ t) %ac f{t)> Eq. ( |T8| ) 
can be expressed as Euler's hypergeometric differential equation 

(t(l - 1)^ + [2 - (A + + A_ + l)t] j t - A + A_) f{t) = , (19) 

with solution f(t) = 2 i 7 i(A + , A_; 2; t), and where the coefficients X± are defined by 

^ f 1 + iac ± iy/c + a 2 c 2 repulsive potential 
\ 1 + iac ± \fc — a 2 c 2 attractive potential 

Thus, the full solution is Xo — t(l — t) iac 2-F1 (A+, A_; 2; t), up to an irrelevant normalization. 

To compute the phase shift 5q, we are interested in the behavior of xo a s x — > 00 (or t — > 1). In 
this limit, we have 7 

*><*> ^ eW " + r ( 2 ( -A + )r ( 2-l) ^ - **" + « • <21) 



where the phase shift is given by 



1 This follows using the identity 

2 F, (A, B; C; t) = ^g^^"^ ~ ^ 2 F 1 (A, B; A + B C + 1; 1 - t) 
T(C)T(A + B-C) n ^c-a-b 



(1 - tf" A ~ u 2 F 1 (C-A,C-B;C-A-B + l;l-t), 



T(A)T(B) 

which is valid for non-integer A + B — C, and also using 2 F% (A, B;C;0) = 1. 
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To the extent that s-wave scattering dominates, we expect ps 4-7T sin 2 £0/ A; 2 to be a useful 
analytic approximation to the full numerical calculation. On the other hand, when m x vjm^ > 1, 
i > partial waves become important and our analytic result is no longer valid. 8 

The existence of s-wave resonances can be inferred from Eq. ( |2~2~| ) by considering the zero 
velocity limit (since s-wave resonances correspond to bound states at zero energy). First, we 



consider the attractive case. Expanding Eq. ([22]) for small a (recall 2a = v/ax), we have 



6 — ► -r2 7 + V(l + \/e)+V'(l- y/c)]ac (23) 

with digamma function ip(z) = T'(z)/T(z) and Euler-Mascheroni constant 7. Thus, as v — > 0, 
the phase shift scales as 5o oc v and ot approaches to a constant. However, this expansion breaks 
down when y/c = n, where n is a positive integer, due to poles in the gamma function. In this 



case, Eq. ( |2~2~| ) gives a maximal phase shift 8 — > ±ir/2 for v — > 0, corresponding to a resonance 
where the cross section is enhanced as ot oc 1/v 2 . In terms of physical parameters, the resonance 
condition is 

°^=n\ ,1 = 1,2,3,... (24) 

As expected, this is the same resonance condition derived for Sommerfeld enhancements ll87l . 
since the same bound state formation is relevant for both scattering and annihilation. We also note 



the appearance of antiresonances (5 = 0), with vanishing s-wave cross section. From Eq. ( |23] ), 
the antiresonance condition is 

ax77lX = r 2 , r« 1.69, 2.75, 3.78, 4.80, 5.81, ... (25) 

where r corresponds to positive roots of the equation 27 + ^(1 + r) + ip(l — r) = 0. On the other 
hand, for a repulsive potential, we have 

8 — > -[27 + ^(1 + iy/c)+^(l- iy/c)]ac. (26) 

^ — >o 

As expected, there is no possibility of resonances, since poles of the gamma function are along the 
real axis only, nor antiresonances, since the quantity in brackets is strictly positive. 

The numerical value of k can be determined a posteriori. In computing the Sommerfeld- 
enhanced annihilation cross section, Ref. [1871 fixed k = ir 2 /6 ~ 1.64 in order to match the 
perturbative result in the Born limit at zero velocity. Applying this prescription to scattering, we 
wish to relate the Born cross section in Eq. (|6]) to our result from the Hulthen potential for v — > 0. 
In the perturbative limit, Eqs. < [23| ) and ([26j) give 5 = ±2£(3)ac 2 , and we have 

Born ^a 2 x m 2 x Hulthen 167Ta|m|C(3) 2 

®T = 4 ' a T = 4—4 • ( 27 ) 

1712 K 1712 



Equating these cross sections gives k = v/2C(3) ~ 1.55. 9 However, since there is no unique exact 



8 Ref. [87 1 generalized this method to £ > by approximating the centrifugal term by a different function allowing 
a solution to the Schrodinger equation. However, the modified centrifugal term alters the long distance behavior of 
the wave function, and the I > phase shifts we obtain by this method do not agree with our numerical calculations. 

9 The small difference in k stems from a difference in matching the Yukawa and Hulthen wavefunctions at r — > 
or r — >• 00. Ref. Il87l obtains k = 7r 2 /6 by equating the wavefunctions at r —> 0, requiring that the integral 
J °° dr' r' V(r') is matched between the Yukawa and Hulthen potentials, using the Lippmann-Schwinger equation. 
Following the same argument, but for r — > 00, one requires J °° dr' r' 2 V(r'), giving our result k — ^/2C(3). 



18 




m^m x matrix 

FIG. 7: Numerical calculation (black solid) and analytic s-wave result (red dashed) for otvt? x as a function 
of m^/mx, for v = 10 km/s, ax = 10~ 3 — 10 _1 , and both attractive and repulsive interactions. The 
analytic approximation breaks down for mxv/m^ > 1, when £ > partial waves become important. 

value for k outside the Born limit, we take simply k = 1.6 which provides an accurate choice 
across a wide parametric range. 

Next, we compare our analytic result for a T with our numerical calculation, shown in Fig. |7} 
Taking a typical dwarf velocity v = 10 km/s, we plot a T m x as a function of m^/mx, calculated 
numerically (black solid) and analytically from the Hulthen potential, with k = 1.6 (red dashed). 
Each panel shows a different coupling, ax = 10 _1 , 10~ 2 , 1CT 3 , for either an attractive (left) 
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ar/mx at v = 10 km/s mx ax 

1 cm 2 /g 210 GeV 2.8 GeV 2.3 x 10~ 2 

10 cm 2 /g 100 GeV 0.67 GeV 1.1 x 10" 2 

TABLE I: Benchmark points for resonant s-wave scattering with da/dQ. = (mxv/2)~ 2 and ot 
16-7r(mxv)~ 2 , consistent with correct DM relic density. 



or repulsive (right) interaction. 10 Our numerical and analytic results agree remarkably well for 
m^/nix > v ~ 3 x 10~ 5 , where scattering is predominantly s-wave, accurately mirroring the 
pattern of resonances and antiresonances within the resonant regime (for the attractive case). This 
agreement provides a highly nontrivial confirmation of our numerical calculation. For m^/mx < 
v, the results diverge as £ > partial waves become more important, as expected. 

Lastly, we provide a series of benchmark parameters for resonant s-wave scattering. This case 
provides a simple and novel velocity-dependence for DM scattering, cr r oc v~ 2 , and it would 
be interesting to incorporate this case within N-body simulations. On-resonance (5 = 7r/2), 
the differential cross section is da/dVL = (mxv/2)~ 2 , giving <jt = I6nm x 2 v~ 2 . In Table I we 
list benchmark parameters that give resonant scattering and also produce the correct DM relic 



density (via p-wave annihilation XX — > cfxp, see Sec. VII). We consider several values of mx 



to fix a T /m x on dwarf scales; the remaining parameters (m^, ax) are determined by resonance 



condition, Eq. ( |24| ) with n = 1 and k w 1.7, and the relic density. For these parameter points, we 
have checked that m x v /m^ < 1 up to cluster scales v ~ 1000 km/s, validating neglect of i > 
partial waves. 



VII. RELIC DENSITY 

In the above discussion, we have taken a x to be a free parameter. In this section, we fix a x for 
a given (m^, mx) through the DM relic density, set by the annihilation process XX — > (jxp. We 
consider here two representative cases where is a vector or scalar field. The annihilation cross 
sections at tree-level are 



/ Ntree 7T«X L m \ ( Ntrcc 3 1XO? x 2 / mj 

m x y m x a 4 m x y m x 

for the vector and scalar mediators. It is clear that DM annihilation to scalar mediators is a p-wave 
process. Since the mediators have masses around 1 — 100 MeV, they will also lead to Sommerfeld 
enhancements for DM annihilation lfi~3l l83l . These enhancements can be important in the early 
Universe for heavy DM. 

The formalism for the symmetric freeze-out with s-wave Sommerfeld enhancements has been 
discussed (84). Here, we expand it to include the p-wave case. 11 The coupled Boltzmann equations 

10 The quantity aTm x is useful to consider since it depends only on m^/mx, after a x and v are fixed, rather than 
m x and separately. Thus, for every point along the curves shown in Fig. [7] one can fix gt / mx to any desired 
value (e.g., 1 cm 2 /g) by taking the appropriate values of mx,m^. 

11 See also [851. 
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for the species X and X can be written as 



dY x ,x [T g* s /^g* 



m pl m x ^^- (a an v) (Y X Y X - Y 2 ), (29) 



dx V 45 x 



where we take the standard definitions 12 x = m x /T and Y xx = n x ^ x /s, with n x X the DM 
number density, s the entropy density, and the equilbrium value of Y x x . In addition m pl ~ 
1.2 x 10 19 GeV is the Planck mass, (a aD v) the thermally- averaged annihilation cross section, and 
g* s and g* are the relativistic degrees of freedom for entropy and energy density, respectively. 

During freeze-out, DM particles have a high velocity and the Sommerfeld enhancement effect 
is negligible. Thus, the freeze-out temperature can be estimated as usual ll86ll 

Xf ~ \n[0.038n(n + l)m p imx(g/y/g^)(To] 

-(n + l) In (In [0.038n(n + l)m pl m x (g/ ^)a }) , (30) 



where g = 2 is the number of degrees of freedom of X and a is given by the relation (a^v) = 
(T x /m x ) n a , where T x is the DM temperature, and n indicates the annihilation type, i.e., n = 
and 1 for s-wave and p-wave annihilation, respectively. 

After freeze-out, F cq becomes insignificant. Neglecting Y cq , we can solve the Boltzmann equa- 
tions analytically as Y xx (x s ) ~ 3.79/(m p \m x J) with 

f Xkd g*s/ \/g* f Xs g*s/ \/g* 

■J = 1 — (a an v)dx+ / — (a an v)dx, (31) 

Jx f x Jx kd x 

where Xkd is the value of x at kinetic decoupling and x s is its value when DM annihilation be- 
comes insignificant and we may stop the integration. Before kinetic decoupling, DM has the 
same temperature as the thermal bath T x = T. After kinetic decoupling at Tkd, the DM velocity 
distribution may be distorted from Maxwell-Boltzmann in scenarios with Sommerf eld-enhanced 
annihilation, since annihilations preferentially deplete the low velocity population. But as shown 
in ||84ll , DM self-interactions mediated by can maintain kinetic equilibrium in the parameter 
region we are interested in, and in this case, we simply take the Maxwell-Boltzmann distribution 
with T x = T 2 /T kA . 

When the DM distribution is thermal with temperature T x , the thermally-averaged cross section 
in the nonrelativistic limit is 



M= / e-W*w (32) 



where v = ^/2T x /m x = ^2/x x . We write the annihilation cross section as a aa v = S(a eLI1 v) tme , 
where (cr an f ) tree is the cross section calculated at the tree-level and S is the enhancement factor. 
Thus, the thermally-averaged annihilation cross section is 



3/2 » 

{ a^v) = / S(a an v) tTee v 2 e- Xxv2 / 4 dv. (33) 



The reader should not be confused with x = axmxr defined in Sec. IV. 
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FIG. 8: The value of ax required to obtain the correct DM relic density as a function of the DM mass mx 
(solid red) for the vector (left) and scalar (right) mediators. We also plot the required ax (dashed blue) if 
the Sommerfeld effect is neglected in the early Universe. We take the DM kinetic decoupling temperature 
Tkd = 1 MeV and the mediator mass = 10 MeV. 

In the cases we consider, the tree-level annihilation cross sections are given by Eq. ( [28] ) and the 
Sommerfeld enhancement factors for s-wave and p-wave annihilations are 

s = * sinh(27rac) ^ = (c - l) 2 + 4(ac) 2 g 

S a cosh(27rac) - cos^^c - (ac) 2 ) ' P l + 4(ac) 2 

respectively, where we have used a = vjlotx and c = 6b/n 2 = 6ax r mx/'K 2 m ( j } [87|. 

In Fig. [8} we show the value of ax which gives rise to the observed relic density for the vector 
(left) and scalar (right) mediators. In the calculation, we have taken the DM kinetic decoupling 
temperature Tkd = 1 MeV and the mediator mass = 10 MeV. The Sommerfeld effect in the 
early Universe can lead to an 0(1) suppression factor on ax for mx > 1 TeV, but is negligible for 
lighter DM. This is because heavier DM requires a larger ax which results in a larger enhancement 
factor on DM annihilation in the early Universe. 

Here, we comment on the dependence of the result shown in Fig. [8] on m^ and T kd . Since a 
large mass hierarchy between m x and is required for DM to have sufficient self-interactions 
to affect structure formation when mx > 1 TeV, the mediator is effectively massless for the 
Sommerfeld enhancement. Thus the result is not sensitive to m^. The value of ax can also depend 
on Tkd- For a small Tkd, DM particles cool down slowly, which suppresses the Sommerfeld effect. 
However, typically, this dependence is very mild because the the DM annihilation rate becomes 
much less than the Hubble expansion rate before the Universe cools to T k d, even if the annihilation 
is enhanced. In our case, we have checked that ax only changes by less than 3% when we set T k d 
to be 1 GeV. It is worth noting, however, that T k d may play an important role in the resonance 
regime. It has been shown that DM can re-couple to the thermal bath after freeze-out in the 
resonance regime, which leads to a negligible relic density [|84l . This chemical re-coupling effect 
only occurs when T k d is high and parameters have to be highly fine-tuned to satisfy the resonance 
condition exactly. With Tk d = 1 MeV, we have checked that chemical re-coupling does not happen 
and DM has the correct relic density in the resonance regime. 

In Fig. |9j we show the allowed range of {mx, m^) with ax fixed by the relic density constraint 
as shown in Fig. [8] For the vector mediator case (left), both attractive and repulsive interactions are 
present, and we take the average of attractive and repulsive cross sections. In the scalar mediator 
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FIG. 9: Parameter space for self-interacting DM as in Fig. [6] with ax fixed to obtain the observed relic 
density via XX — > <p<p annihilation at freeze-out. The left (right) panel shows the vector (scalar) mediator 
case whwere annihilation is s-wave (p-wave). Crosses show benchmark points in Table VI The lines and 
colored regions are as in Fig. [61 



case (right), DM self-interactions are purely attractive. It is clear that the allowed region for 
solving the small scale anomalies is still broad even after we impose the relic density constraint 
on ax- 



VIII. OBSERVATIONAL TESTS 



Self-interacting DM has distinct signatures in direct detection experiments because self- 
interactions thermalize the DM velocity distribution [|92l . In this section, we discuss signatures 
of self-interacting DM in indirect detection observations, when DM in halos self-annihilates. As 
we have shown, the existence of a light mediator is essential for generating a large enough self- 
scattering cross section. The same mediator can also lead to Sommerfeld enhancements for DM 
annihilation in halos if DM is symmetric. Since the enhancement effect increases as the DM 
velocity decreases, we expect DM particles in dwarf galaxies to have a larger self-annihilation 
cross section than those in the Milky Way or clusters. This scale-dependent feature of the DM 
annihilation cross section can be potentially determined by studying signal fluxes from different 
astrophysical objects. 

to 



Here, we take a few examples from the self-interacting DM models given in Section VI 



show Sommerfeld enhancements for DM annihilation in halos. We consider the case where DM 
particles annihilate to SM states in DM halos with s-wave processes. 13 To illustrate the point in a 



13 A familiar example is usual symmetric DM. Asymmetric DM can also generate annihilation signals if DM-anti-DM 
oscillations occur in the late epoch I8& M9T1 . 
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FIG. 10: The thermally-averaged s-wave Sommerfeld enhancement factor (S) and transfer cross section 
(cry) as a function of mx/m^ for ax = 10 -2 (left) and ax = IO -3 (right) with vq = 10 km/s (blue), 
200 km/s (red) and 1000 km/s (green), corresponding to the most probable DM velocities on dwarf ("dw"), 
Milky Way ("MW"), and cluster ("cl") scales. One can see the correlation between the enhancement in the 
annihilation cross section and the scattering cross section due to the s-wave resonance. 



rather model-independent way, we take the assumption that DM has the correct rel ic de nsity and 
do not demand XX — >■ 



VII 



We have 



to set the correct relic density as discussed in Section 
checked that our result does not change qualitatively if we demand the relic density set through 

.V .V -> oo. 

For s-wave annihilation, the relative annihilation rates on different scales are determined 
by Sommerfeld enhancements folded together with DM distributions. Of course, DM self- 
interactions will also alter the density profiles in the center of the DM halos, changing the an- 
nihilation rates. Rather than folding the DM distribution in to extract the total rate, we focus on 
the effect of the Sommerfeld enhancement alone on the annihilation cross section. We calculate 
the thermally-averaged Sommerfeld enhancement factor as 



(S) 



(27TO 2 ) 3 / 2 



(35) 



where S s is the s-wave Sommerfeld enhancement factor given in Eq. ( |34| ). 
In the top two panels of Fig 



10 we plot the thermally-averaged Sommerfeld enhancement 



factor for DM annihilation as a function of mx/m^ for different ax and DM velocities. The 
upper two panels in Fig. [10] are complementary to those in Fig. [6] and to the lower two panels of 
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Fig. [10} which show the preferred parameter space for the self-scattering cross section to solve 
the small scale structure problem. One can see the correlation between the enhancement in the 
annihilation cross section and the scattering cross section due to s-wave resonances. It is also 
clear that, similar to the DM self- scattering case, there are three distinguishable regions for the 
Sommerfeld enhancement factor depending on mj/m^. If the mediator and DM masses are 
comparable, it is in the Born regime where (S) is negligible on all scales. On the other hand, if the 
DM mass is much larger than the mediator mass, the enhancement factor becomes independent 
of mx/rritf, which corresponds to the Coulomb limit. In this limit, the enhancement factor is 
essentially given by S ~ nax/v. In the third region where mx/m^ ~ n 2 n 2 /6ax with n = 
1, 2, 3..., DM annihilation can be enhanced resonantly. On resonance, the enhancement factor is 
S ~ ir 2 ax'm ( f > / (6mxv 2 ) which is very sensitive to the DM velocity. We emphasize that the s-wave 
resonant condition for the Sommerfeld enhancement of DM annihilation exactly corresponds to 
DM s-wave resonant self- scattering. 

As shown in Fig.|6j most of the parameter space preferred for solving the small scale structure 
problem is in the resonant and classical regions. In these regions, constraints on the DM self- 
interacting cross section from DM halo shapes and the Bullet Cluster are elegantly evaded by the 
velocity-dependence of the self-scattering cross section. Interestingly, in the same regions, the 
Sommerfeld enhancements for the DM annihilation cross section differs significantly on different 
scales. In the resonant (classical) region, (S) in dwarves can be a factor 0(100) (O(10)) larger 
than that in the Milky Way. Therefore in many cases self-interacting DM predicts Sommerfeld 
enhancements for DM annihilation. If indirect detection signals are observed and annihilation 
cross sections are measured on different scales, it will give us a strong hint for self-interacting DM 
and help us further narrow down the parameter space. 

IX. CONCLUSIONS 

We have examined DM self-interactions via a Yukawa potential with a massive dark force. Over 
much of the parameter space, the Born (ajini < m^) and classical (mxv > m^) analytic for- 
mulae break down, and quantum resonant structures, many with non-trivial velocity or angular de- 
pendences, arise. We devised a method that allowed us to efficiently explore the strongly-coupled 
regime of parameter space. We examined in detail the structure of this regime, and matched our 
results onto the known classical formula, verifying for the first time that analytical result. We were 
also able to derive an analytic formula for our results for the case of a strongly-coupled s-wave 
resonance. We also extracted the angular dependence of our results in the quantum and classi- 
cal regimes, adding another dimension for study to the dynamics of DM self- scattering which is 
particularly important when the mediator is light. 

Our results have implications for the future study of DM self-interactions. Theoretical study 
and simulations of DM self-interactions have focused on simple analytic solutions for the scat- 
tering cross section, with constant or classical velocity (and no angular) dependence. New sim- 
ulations are in progress which will better account for baryonic effects on DM structure, while 
simultaneously integrating DM self-interactions ||93l . It will be important to simulate a broader 
class of DM self-interaction models by including strongly-coupled and resonant effects in the 
simulations. Angular dependence should also be modeled, though including the general angular 
dependence in the strongly-coupled regime can be difficult. However, we found a few cases where 
the scattering cross section has the desired velocity-dependence while the angular dependence is 
rather simple. In the case of s-wave resonant scattering, the scattering cross section scales as v ~ 2 
and is also isotropic. In the strongly-coupled classical regime, we have numerically confirmed 
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that isotropic assumption for scattering on dwarf galaxy scales which has been taken in the recent 
simulation In addition, the Rutherford formula is available in the massless mediator limit. We 
have devised benchmarks which may be utilized in simulations. 

In addition, our results allow the correlation of DM self- scattering with annihilation, having 
implications for indirect detection experiments. Sommerfeld enhancements for DM annihilation 
directly correspond to velocity dependent self-interacting DM. Conversely, the absence of Som- 
merfeld enhancements imply a velocity-independent DM self- scattering cross section, so that if 
cores form in dwarves they also form in clusters. 

Clearly DM self-interactions provide an avenue for exploration with rich consequences for DM 
structure in our Universe. While the nature of the DM may first be revealed through its interactions 
with ordinary matter, to date everything we have learned about DM has been gleaned through the 
formation of structure. DM self-interactions can change this structure in complex ways, so that as 
we learn more about it, we may also uncover evidence for the particle physics nature of DM. 
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Appendix A: Compendium of analytic results and benchmark points 

We summarize analytic results for self-interacting DM scattering through a Yukawa potential. 
The relevant parameters are the DM mass mx, the dark force mediator mass and coupling qj, 
and the relative velocity v. The transfer cross section ot = J dQ(l — cos6)da / d£t provides a 
useful proxy for comparing specific particle physics models to N-body simulation results. We also 
give do /dVL, which is a required particle physics input for simulations. 

In the Born limit (axmx/m^ <C 1), the cross section can be computed perturbatively in ax- 
The differential cross section is do/dVL = a x m x /(m^ + m x v 2 (l — cos#)/2) 2 , giving 



On 



.B» _ 8™| / (1 + m2 2 2) _ _^_ ) (A1) 



m x v 4 V ml + m x 



for both attractive and repulsive potentials ETTl . 

non-perturbative effects become important outside the Born regime (axmx/m^ > 1). Results 
have been obtained in the classical limit (mxv / ^> 1), giving for an attractive potential [|2Tl[80l 



rnh 1 



a 



-■las J 8^2/ ^ + L5/3 1.65) 1Q -1 < p < 1Q 3 ^ 



^(ln/^ + l-im- 1 /?) 2 $> 10 3 
and for a repulsive potential [|20l [811 



Mln2/3-mln2/3) 2 fi>l ^ 
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where (3 = 2axm ( j } / {mxv 2 ). We find that da/dVt ps crr/(47r) (i.e., approximately constant) 
for {3 < 1, but approaches the Rutherford scattering formula da/dVt ps a x /(m^t> 4 sin 4 0/2) for 

> 1. 

Outside the classical regime (mxv/m^ < 1), the cross section is largely dominated by s-wave 
scattering. We have obtained a new exact non-perturbative result for <j T for the Hulthen potential, 
which provides an excellent approximation for the true Yukawa potential. Our result is: 

ff Hulth*n = J67T ^ (M) 

where the £ = phase shift is given in terms of the T-function by 




£o = arg , A± = < v r- (A5) 



^ ± /axmx_^|4 attractive 



l + *mxv ±i I ffla + repulsive 



and k p^ 1.6 is a dimensionless number. The differential cross section is da /dVL = <tt/(47t). This 
formula takes into account non-perturbative effects associated with s-wave scattering, and covers 
a complementary parameter region to the classical and Born formulae. 
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